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By means of numerical simulations and epidemic analysis, the transition point of the stochastic, 
asynchronous Susceptible-Infected-Recovered (SIR) model on a square lattice is found to be co — 
0.1765005(10), where c is the probability a chosen infected site spontaneously recovers rather than 
tries to infect one neighbor. This point corresponds to an infection/recovery rate of Xc ~ {1 — 
Co)/ Co = 4.66571(3) and a net transmissibility of (1 - Co)/(l + 3co) = 0.538410(2), which falls 
between the rigorous bounds of the site and bond thresholds. The critical behavior of the model 
is consistent with the 2-d percolation universaUty class, but local growth probabilities differ from 
those of dynamic percolation cluster growth, as is demonstrated explicitly. 



I. INTRODUCTION 

The Susceptible-Infected-Recovered (SIR) model is a 
fundamental system in epidemiological modeling [IHH] 
and has been studied extensively within the context of 
non-equilibrium phase transitions and critical phenom- 
ena (i.e., [nHlSlIU])- The model was developed to de- 
scribe the propagation of an epidemic that occurs during 
a period of time much smaller than the lifetime of indi- 
viduals of a given population. It is assumed that the pop- 
ulation can be divided into three categories: Susceptible 
(S), Infected (I) and Recovered (R) individuals. Suscep- 
tible individuals become infected at a given rate through 
contact with infected individuals. Infected individuals 
recover with a given rate and become immune and recov- 
ered. The model is capable of showing a threshold of the 
epidemic spreading as one increases the infection rate. 

The SIR process has been studied using different ap- 
proaches and contexts. Originally it was defined in 1927 
by Kermack and McKendrick [1 as a deterministic pro- 
cess by means of a set of ordinary differential equations; 
they showed that epidemics disappear before all the sus- 
ceptible individuals contract the disease. Afterwards the 
model was given a stochastic description by means of 
birth and death processes [5], and later Grassberger [5] 
introduced a cellular automaton implementation, that is, 
a synchronous-update Markovian process on a lattice. In 
this paper we consider a stochastic, asynchronous-update 
lattice version of the SIR model fT5^ , in which lattice sites 
are updated one at a time. This model is a special case 
of the predator-prey stochastic lattice-gas model intro- 
duced by Satulovky and Tome [25l [26] , and also consid- 
ered by Antal et al. [T^l [27]. It is also a special case 
of the Susceptible- Infected- Removed-Susceptible (SIRS) 
stochastic lattice gas model [19] . For the synchronous or 
asynchronous versions of the SIR process, one observes 
that, as the model parameters are varied, a phase transi- 



tion takes place. This is a continuous phase transition be- 
tween two distinct regimes: one in which the population 
remains susceptible (non-spreading regime) and another 
in which the epidemic spreads over the lattice (spread- 
ing regime) , where a significant portion of the population 
becomes infected and eventually immune. At the transi- 
tion the system becomes critical and corresponds to the 
epidemic threshold. Cardy and Grassberger JU)^ argued 
that the transition found in the SIR cellular automaton is 
of the percolation universality class. This has been con- 
firmed in various ways by numerous studies in stochastic 
models with synchronous as well as asynchronous update 

[Iii[ia[i5i[ii[i8i[i9]. 

In recent years there has been a great deal of interest in 
the SIR model on networks and other systems (i.e., [131 - 
[TOl [HHH [2HH3S]). In this paper we consider the SIR 
model on the square lattice, which has wide applications 
to many physical problems such as the spread of disease 
in plants, and which has been studied only to a limited 
extent (S] [TH [THl [3S] ■ One question that has come up in 
the context of mean field and network studies is the ef- 
fect of heterogeneity in the infectious period [23 [301 [32] • 
When the period of infection is fixed, there is a direct 
mapping of the model to bond percolation on the lattice, 
and thus the transition point can be determined exactly 
[21 [m [Ml [SZ]- However, when the infection period is 
heterogeneous, such as the exponentially distributed in- 
fection period inherent to the original SIR model, there 
is no exact solution. In this paper, we carry out a careful 
numerical study of the transition point for the exponen- 
tially distributed case, and then investigate the correla- 
tions in the transmission of the infection to the nearest 
neighbors. We work out the correlations explicitly, and 
find that at the critical point there is a higher probability 
to infect more rather than fewer of the neighbors, unlike 
the case of bond percolation where the distribution of 
infected neighbors is simply binomial with p — 1/2. 

The rest of the paper is organized as follows: In section 
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\n\ we define the asynchronous SIR model in terms of the 
master equation. In section III we describe the Monte- 



Carlo algorithm that we use, and in section [IV| we carry 
out a numerical analysis of the cluster size distribution of 
recovered individuals to precisely determine the critical 
point of the model. In section |V] we show in detail how 
the asynchronous SIR model differs from percolation on 
a local scale. In section IVll we discuss the relation of our 
results to some other population biology models, and in 
section [VII| we give our conclusions. 



II. THE SIR LATTICE MODEL 

To model the dynamics of an epidemic with immu- 
nization we consider a stochastic lattice-gas model with 
asynchronous dynamics. The lattice plays the role of the 
spatial region occupied by the individuals and the lattice 
sites are the possible locations for the individuals. Each 
site can be occupied by just one individual that can be 
either a susceptible, an infected, or a recovered individ- 
ual, called, respectively, an S site, I site and R site. At 
each time step a site is randomly chosen and the follow- 
ing rules are applied: (i) If the chosen site is in state S 
or R it remains unchanged, (ii) If it is in state I then 
(a) with probability c the chosen site becomes R and (b) 
with the complementary probability b = 1 — c a neighbor- 
ing site is chosen at random. If the chosen neighboring 
site is in state S it becomes I; otherwise it remains un- 
changed. Notice that a site in state R remains forever 
in this state so that the allowed transitions are S— >I— >R. 
From this set of dynamic rules it follows that the state 
of the system will change only when the chosen site is 
in state I, a feature that will be used to speed up the 
simulation as explained in section |III[ This algorithm is 
equivalent to making S— t-I with probability 6ni/4, where 
ni is the number of nearest-neighbor I sites, and I— )-R 
with probability c. 

The system evolves in time according to a master equa- 
tion. To each site i of a two-dimensional lattice, we as- 
sociate a stochastic variable rji that takes the values 0, 
1, or 2 according to whether the site i is occupied by 
an R or an S or an I individual, respectively. A micro- 
scopic configuration of the entire system is denoted by 
the stochastic vector rj — {rji, . . . , 77^, . . . , ry^v) where TV is 
the total number of sites. Because the possible transi- 
tions are the cyclic ones (1 — 2 — 0), it is convenient 
to define the state ry* obtained by a cyclic permutation 
of the state of site i, that is, 77* = (ryi, . . . , 77^, . . . , iji^) 
where 77- is 1, 2, or according to whether rji is 0, 1, 
or 2, respectively. The master equation for the probabil- 
ity distribution P{ri, t) associated with the microscopic 
configuration rj at time t, is given by 



dt 



(1) 



where the summation on i extends over the nearest- 
neighbor sites of site j and ^rj denotes the state obtained 
from 77 by an anticyclic permutation of state (0— ^2— 
The quantity wfj is the transition rate associated with 
the infection process, given by 



/3 



(5(77„1)(5(77„2) , 



(2) 



where z is the number of nearest-neighboring sites of site 
j (the lattice coordination number) and S{x,y) is the 
Kronecker delta; and is the transition rate for the 
recovery process, given by 



(3) 



Two external parameters are associated to these pro- 
cesses: the infection rate /3 and the recovery rate 7. 
The probabilities b and c are related to these rates by 

b = (i/{l3 + 7) and c = 7/(/3 + 7) so that 



6 + c= 1 



(4) 



as it should. 

From the master equation ([T]) we can derive the time 
evolution equations for the densities of recovered po, sus- 
ceptible pi, and infected p2- The connection of the 
present stochastic lattice model with the approach de- 
veloped by Kermack and McKendrick |lj can be revealed 
by using a simple mean-field approximation. Within this 
approximation [19' the following set of ordinary differen- 
tial equations for the densities can be derived 



dt 



pi = -Ppip2 , 



dt 



dt 



(5) 



(6) 



(7) 



These equations are essentially the equations introduced 
by Kermack and McKendrick [T] in their deterministic 
approach for the spreading of a disease with immuniza- 
tion. 

To analyze an epidemic, one can begin from an initial 
condition at time t = where all the individuals are sus- 
ceptible, with the exception of a very small number of 
infected individuals; that is, pi = 1 — p*, p2 = P* ^ 1 
and po = 0. Using this initial condition and Eqs. ([5][7|, 
one finds that the system evolves in time and reaches two 
types of states: one where the epidemic spreads, that is, 
po ^ 0, pi ^ 0, p2 = 0, when t — >■ 00, which occurs 
for sufficiently large values of the infection probability b, 
and another where the epidemic does not spread, that is, 
pi = 1 when i — > cxD, which occurs for small values of 
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b. As one varies the parameter b, there is a continuous 
phase transition at 5 = 1/2. In stochastic lattice models 
one expects similar behavior but a distinct value for the 
critical parameters. In that case one starts with a lat- 
tice full of susceptible individuals with the exception of a 
single infected individual, and studies the epidemics that 
ensue. 



III. SIMULATION ALGORITHM 

The asynchronous SIR model may be simulated by a 
kinetic Monte-Carlo process by following procedure: 

• Pick an I site randomly from a list of all I sites. 

• Generate a random number X in (0, 1). li X < c 
then let I — >■ R. 

• Otherwise (if X > c), pick one nearest neighbor of 
the I site randomly. If the neighbor is S, then make 
it I and add to list of I sites. 

• Repeat as long as there are available I sites. 

When we remove an I site from the list in the first step, 
we swap the empty location with the site at the top of 
the list and decrease the list length by one in order to 
keep the list compact. New I sites in the third step are 
added to the top of the list. 

The Monte-Carlo time t is determined by increment- 
ing t by where A^i is the current number of I sites 
on the list, each time an I is picked from the list. This 
reflects an effective increase of one Monte-Carlo step if 
every site on the lattice were chosen once, on the aver- 
age. However, as explained below, we don't actually keep 
track of time in our simulations, but instead monitor just 
the size (number of R sites) of the clusters. 

This model is closely related to standard percolation 
growth, in which an active growth site spreads to an un- 
occupied nearest neighbor with the given bond probabil- 
ity. However, in the percolation described below, 
the growth site becomes inactive after the four nearest 
neighbors are checked (for both asynchronous and syn- 
chronous updating). In the SIR model, the I site re- 
mains active an exponentially distributed length of time, 
and can thus attempt to infect neighboring sites multiple 
times. That behavior leads to a kind of correlation in the 
spreading in the SIR model, as we shall see. 

In order to study properties of this model in detail, it is 
necessary to know the transition point to high accuracy. 
Recent work [TS] has shown that cq = 0.1765(5), where 
the number in parentheses represents the error in the last 
digit. In this, paper, we find cq to about 200 times the 
accuracy, and indeed, in so doing, we find nearly complete 
agreement with the previous central value. 



IV. NUMERICAL RESULTS 

To find the transition point, we consider the statistics 
of the cluster size distribution, and determine cq as the 
point where the distribution follows a power law. Ac- 
cording to standard percolation theory, the probability 
Ps that a point belongs to a cluster containing s sites is 
given by sug , where Ug is the number of clusters of size s 
(per site) on the lattice. At criticality, Ug ^ s~^, so that 
Ps ~ s^~'^. Integrating from s to oo we find the proba- 
bility P>s that a point belongs to a cluster of size greater 
or equal to s. At criticality, P>s s^^^i a-nd within the 
scaling region one expects (SB] 

P>s - s^-^F{es'') « s^-^{A + Bes") , (8) 

where e ~ c — cq and for the last term we made a Taylor- 
series expansion of the scaling function F for small val- 
ues of its argument [39] . In 2-d percolation, we have 
r = 187/91 and a = 36/91. While there is no "n^" 
for dynamic percolation (because we don't have a lattice 
fully populated with clusters), we expect P>s for dynamic 
percolation and the SIR model at criticality to show the 
same power-law behavior as static percolation. 

To calculate P>s, we ran simulations with the lattice 
covered entirely by S sites, except for a single I site in 
the center, using the algorithm above. The lattice size 
was 16384 x 16384. The size s of the cluster was charac- 
terized by the number of R sites. When s went beyond 
a cut-off of 2^1 = 2097152, the growth was stopped. The 
cluster sizes were binned to make a histogram. Various 
runs were made at each value of c as listed in the caption 
of Fig. [T] For a pair of values of c which bracketed the 
apparent transition point, we generated about 10^ clus- 
ters, requiring several weeks of computer time each. The 
random number generator we used was R(9689) of [10] . 

While the time t might perhaps be more of a natural 
variable to consider for dynamic percolation, we chose to 
consider the survival probability as a function of s. One 
advantage of using s instead of t is that the percolation 
exponents for P>s are known exactly in two dimensions, 
while those for P>t ~ t(2-T)£'/<i„.„ ^ ^-0.09213... 
icality) are related to dmin which is known only approx- 
imately: d,„i„ « 1.1306(3) [H]. (Here, D = 91/48 is the 
fractal dimension.) Also, it is somewhat easier to keep 
track of the dependence upon a discrete variable s rather 
than a real variable t. In Fig. [l] we plot s'^^^P>s vs. s'^ 
for various c, using the 2-d percolation values of r and a. 
According to ([8|, this plot should be a linear function, 
with a slope proportional to c — cq, and slope equal to 
zero for c — cq. Indeed, this linear behavior is followed 
very well, except for smaller s where ([s]) is not valid due 
to finite-size effects, in Fig. [2j we plot the slopes of the 
linear parts of the three central values as a function of 
c. The intercept where the slope is zero gives the critical 
point: 

Co = 0.1765005(10) (9) 
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FIG. 1. (Color online) Plot of s^-2p>, vs. s" for the 
asynchronous SIR model, using dynamic percolation values 
r = 187/91 and a = 36/91, for (top to bottom for large s): 
c = 0.176490, 0.176495, 0,176500, 0.176505, and 0.176510. 
The number of samples are, respectively, 335000, 10500000, 
440000, 11000000, 360000. 



where the error is based upon the statistical error of the 
data. The linear behavior of the curves in Fig. [T] for 
larger s shows that the critical behavior is consistent with 
percolation scahng. In Fig. [3j we plot In P>s vs. Ins at 
c = 0.1765, and find a slope —0.05524 consistent with 
the percolation prediction 2 — r = —5/91 ~ 0.054945. 

We mention that we also confirmed that cq ~ 0.1765 is 
the transition point for the alternate simulation method 
(as used in [T^]) in which all sites are sampled, and I— >R 
with probability c and S— >I with probability ni{l — c)/4. 
This simulation runs more slowly because of the time 
spent testing S sites with no I nearest-neighbors (al- 
though it can be speeded up by using a list of eligible 
S sites). The program can also be speeded up somewhat 
by raising both probabilities the maximum amount: I-^R 
with probability c' = c/(l — c) and S— >I with probability 
ni/4. We also confirmed by simulation that this proce- 
dure yields the correct critical point c' = co/(l — cq) ~ 
0.21433. 

In terms of the probability of infection (or the 
net transmissibility) our result implies by (12), 
prob(infection) = 0.538410(2), which falls within the rig- 
orous lower and upper bounds [91I36] of pi''™''^ 1/2 
0.592746. 
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V. RELATION TO PERCOLATION GROWTH 

While the asynchronous SIR model is clearly in the 
dynamic percolation universality class, it is not strictly 
identical to percolation. In this section we show explicitly 
how the two differ locally by considering the probabili- 
ties of an infection spreading from a single I site. This 
difference between the SIR model and bond percolation 



FIG. 2. Slope of the linear part of the three central curves of 
Fig. [l] measured for s"^ > 125. The intercept of the line on 
the horizontal axis gives an estimate for cq. 



was noted by Kuulasmaa [35] and by Grassberger |5] , and 
more recently discussed in general by Kenah and Robins 
[33] . Here we carry out a brief analysis directly related to 
the computer algorithm we developed, and give explicit 
numerical results for the infectivity or transmissibility at 
the critical point. 

First, for comparison, we formulate epidemic bond per- 
colation growth in the SIR language. As in the algorithm 
in section [nTj we start with a system of all S, with one I 
site at the center. The algorithm we follow for percola- 
tion is: 

• Pick an I site randomly from the list of I sites. Set 
this I to R and remove it from the list. 

• Consider the I site's four nearest-neighbor sites, 
and for each do the following: 

• If the neighbor is an S site, then generate a random 
number X in (0, 1) 

9 li X < p then let the S become I, and add to list; 
otherwise, do nothing. 

• Repeat as long as there are available I sites. 

This is equivalent to bond percolation because each 
possible bond is considered only once and with a fixed 
probability p. Note that the algorithm does not generate 
all the bonds of a cluster, as no internal bonds are consid- 
ered, but it finds all the "wetted sites" with the correct 
probability. The bonds that are occupied form a mini- 
mum spanning tree on the cluster. A similar algorithm 
was given by Grassberger [5]. 

If a given I site has m nearest-neighbor S sites (to — 

0, ...,4), then the probability P^™' that exactly k of 
them become infected as a result of the central I site 
for percolation is just 



i—k 



(10) 
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At the threshold Pc = 
havePg^*^ = 1/16, Pi^'^ 

5(4) 



1/2, for m = 
= 4/16, 



4, for example, we 
= 6/16, P|^^ = 4/16 



and Pi = 1/16, with an average occupancy of 2. 

Next we calculate P^'""* for the SIR model. First, we 
derive the probability of spreading for "transmissivity" 
for the SIR model. Consider an I site surrounded by at 
least one S site. We desire the probability that a given 
one of those S sites will become infected by the I site. 
The I site will remain infectious for an exponentially dis- 
tributed number of trials. The probability p{n) that it 
remains infectious for n trials is simply 



p{n) = (1 - cf 



0,1,2,.. 



(11) 



because c is the probability that it recovers in a given 
trial. This implies that, on the average, an I site will 
be considered X)^o ~ (1 " c)/c times before it 
recovers. 
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FIG. 3. Plot of lnP>ji vs. Ins for the SIR model, simulated 
on a 16384 x 16374lattice, with a cutoff of s = 2097152. 
The equation represents a linear fit of the points for s > 100, 
where x represents Ins and y represents lnP>s. 

For each trial that occurs on the given I site, we check 
one of its four neighbors randomly to see if any is S, not 
remembering if we have already checked that neighbor 
before. The probability that a given nearest neighbor 
S site is chosen in at least one of the n trials is given 
by 1 - (3/4)". Multiplying by p{n) from Eq. ([iT]) and 



summing, we find the net probability a given nearest- 
neighbor S site becomes infected: 



Prob (infected) = JZ^^ " 



n=0 
1 - C 

l + 3c 



1 



(12) 



If we set this probability equal to 1/2 as in standard 
bond percolation, we would find c = 1/5, which is above 
the observed value Cq = 0.1765. Thus, finding where 
the effective bond probability equals 1/2 does not give 



the correct threshold. This is because that bond is not 
occupied independently, but is correlated with its neigh- 
bors, and so the effective bond probability (probability of 
spreading to a neighbor) is not 1/2. Putting c = 0.1765 
into this formula, we find Prob (infected) = 0.5384, which 
is greater than in standard bond percolation. 

In another limit, the SIR model becomes identical to 
site percolation on the square lattice [9j. Here, one as- 
sumes that if an I site does not recover in its first trial, 
all neighboring S sites become infected. The thresh- 
old corresponds to 1 — p(0) = 0.592746 — p^c^^'^\ or 
c — \ — — 0.407254, where pt^'^'^^ is the site perco- 

lation threshold on the square lattice p2Vl44j . and have 
used p{{)) from Eq. ( 11 ). This also implies Prob (infected) 
= 0.592746. Thus, we see that the net rate of infection of 
the SIR model (0.5384) falls between that of the (uncor- 
related) bond percolation and the highly correlated site 
percolation values, as mentioned earlier. 

Another comparison to make is to the SIR model 
defined on the Bethe lattice (Cayley tree). Say that 
we have a Bethe lattice with coordination number z. 
Then an extension of equation ( |12[ ) to arbitrary z gives 
Prob(infected) = h/[z — [z — l)6[~^Setting this equal to 
l/(z — 1), because an I site has (z ~ 1) nearest-neighbor 
S sites and the critical point is when, on the average, one 
nearest-neighbor site is infected, we find [53] 



b = 



2(z-l) 



(13) 



Thus, for z = 3, the transition is at c = 1/4, and for 
z = 4, the transition is at c = 1/3. The latter is much 
greater than the critical cq we found for the square lattice. 
This is expected, because for the Bethe lattice the net 
bond probability should be less than for a regular lattice 
with the same z. 

Now, we turn to the probabilities P^™"* for the SIR 
model. From a recursion relation analysis, and using the 
generating function method, we find the coefficients for 
each m = 1,2,3,4. The derivation is given in t he Ap- 
pendix. The resulting P^"*'' calculated from Eqs 



in the Appendix at c = 0.1765 are given in Table Jj I'he 
distribution is seen to be much different from a binomial 



28 



29) 



distribution of standard percolation Eq. ( 10 ), with an in- 



crease in the probability with increasing k for a given to. 
The probabilities add to 1, and the mean number of in- 
fected neighbors is 6/(4 — 36)to = 0.538411to, consistent 
with Eq. ( [T2| ). 

For to = 4, the results give the outgoing probabil- 
ities in the construction of Kuulasmaa [36^, in which 
a quenched double-directed lattice showing all possible 
spreading probabilities is considered. In this construc- 
tion, one does not consider whether the neighbor is S or 
not, but draws outgoing directed bonds (with the proper 
probabilities) to all neighboring sites. Then an epidemic 
is identified by following all directed paths emanating 
from a given seed. 

Thus, we see in general that in the SIR model at crit- 
icality there is an enhancement in the number of neigh- 
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TABLE I. P^™' for the SIR model at c = 0.176500. For 
comparison, the last row (4B) shows the binomial distribution 
of bond percolation at pc = 1/2, which is also at a critical 



point but with much different P^: 
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FIG. 4. (Color online) Example of a completed bond percola- 
tion cluster of 51456 wetted sites, generated at the threshold 
Pc = 1/2. The red (dark) dot marks the place where the 
cluster growth began. 




FIG. 5. (Color online) Example of a completed SIR cluster of 
51034 recovered sites, generated at the threshold co = 0.1765. 
The red (dark) dot shows the location of origin of the infec- 
tion. 



bors that are infected, compared to the case for dynamic 
percolation. Because of the correlations in the infected 
neighbors in the SIR model, it is necessary for the net 
neighbor infection probability to be higher than the ran- 
dom percolation value of 1/2. 

We note that for c — 1/5, the value obtained by setting 
the probability of spreading ( 12 ) equal to the bond perco- 



lation value 1/2, -P^™'' is simply a uniformly distribution 

= 1/(to+ 1). This behavior is in sharp contrast to 
( 10 1 for p = 1/2 or indeed any value of p. 



For the z — A Bethe lattice at its critical point c — 1/3, 
the values of ^ are 5/15, 4/15, 3/15, 2/15 and 1/15 
for fc = 0, 1, 2, 3, and 4, respectively — trending in the 
opposite direction as for the square lattice. 

In Figs. [3] and [5] we show actual pictures of critical clus- 
ters containing about 50,000 sites each from the percola- 
tion and SIR models, respectively. On a local scale the 
SIR clusters appear slightly denser but otherwise there is 
no apparent difference between the two critical clusters. 



VI. RELATION TO OTHER POPULATION 
BIOLOGY MODELS 

The SIR model can be considered as a particular case of 
other population biology models such as the Susceptible- 
Infected-Recovered-Susceptible (SIRS) model [HI IH] 
and the Predator-Prey model [ig. The SIRS model de- 
scribes an epidemic process without permanent immu- 
nization and is defined by the following three processes: 
S— )-I, I— )-R, and R— )-S. The first two processes, S— >-I and 
I— T^R, are the same as those for the SIR model, as de- 
scribed in section 2, and occur with rates /? and 7, re- 
spectively. The third process R— ;>S is spontaneous and 
occurs at rate a. The Predator-Prey model, in the epi- 
demic language, is similar to the SIRS model except that 
the third process, R— >-S occurs with rate ans/4, where ns 
is the number of nearest-neighbor S sites. In the epidemic 
language the following correspondence is used: prey as S, 
predator as I, and an empty site as R. 

The SIRS and the Predator-Prey models exhibit non- 
equilibrium phase transitions between an absorbing sus- 
ceptible phase and active phase where the individuals are 
continuously being infected. Their critical behavior be- 
longs to the universality class of directed percolation for 
a ^ 0. When a = one recovers |TH1 HH], from these 
two models the SIR model, which belongs to the univer- 
sality class of dynamic isotropic percolation, as we have 
confirmed. 

In the opposite regime, namely, when a is large enough 
compared to /3 and 7, both the SIRS and the Predator- 
Prey models map |Tni SS] into the contact process (CP) 
[iB] with a creation rate A = (3/j. The CP can be identi- 
fied as the Susceptible-Infected-Susceptible (SIS) model 
of epidemiological modeling [8j . The SIS model describes 
the dynamics of infection with no immunity and has two 
processes: S— >I with a rate /3ni/4, and I— t^S with rate 7. 
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Changing the time scale, these two processes can be de- 
scribed by the infection rate A = /3/7 and recovery rate 
equal to 1. A site occupied by a particle in the CP model 
corresponds to an I site in the SIS model and an empty 
site to an S site. 

The SIR and the SIS models correspond then to two 
extremal behaviors of the SIRS model with respect to 
immunity. Let us consider the periods of time spent by 
an individual in the R state. These periods of time are 
distributed around a mean value T, which is proportional 
to the inverse of the rate a. In the SIRS model T is finite 
which means that individuals have a partial immuniza- 
tion. The SIR model can be understood as the SIRS 
model in which T — > oo, meaning that an R individ- 
ual has a lifelong immunity. The SIS model, on the other 
hand, can be regarded as the SIRS model in which T — 
implying that an individual has no immunity or equiv- 
alently that an infected individual becomes susceptible 
without passing by the R state. 

Now we would like to give a comparison between the 
critical parameters of the SIR and SIS models defined 
on a square lattice. For this purpose, it is convenient to 
use the parametrization b — j3/{a -f /3 + 7), c — ^/{a + 
/3 -|- 7) and a = a/{a + (5 + so that 6, c and a are 
interpreted as the probabilities of infection, recovery and 
re-infection, respectively. In this formulation there are 
just two independent parameters, as a -\- h + c = 1. 

For small values of /3 and 7 compared to a, the crit- 
ical line of the SIRS model can be obtained by means 
of the mapping into the SIS or CP [121 US] . In two di- 
mensions the critical value for the creation rate for the 
CP is Ac = 1.64874(4) [4H49]. Because A = /3/7 and 
b/c — /3/7, it follows that the critical line is given by 
b/c ~ Xc ~ 1.64874(4) as a ^ 1. If one extrapolates this 
line toa = 0, 6-fc=l, which corresponds to the SIR 
model (that is, suppressing the parameter a), one gets 
(1 - c)/c = Xc or c = l/(Ac -f 1) ~ 0.37754 for the net 
probability that I— >S rather than tries to infect one neigh- 
bor. This value should be compared to the critical value 
for the SIR process cq ~ 0.1765. Alternatively, we may 
compare the infection probability b = 1 — c « 0.62247 
with the critical infection probability for the SIR model 
60 = 1 — Co ~ 0.8235. It follows that epidemic spreading 
for the SIR model occurs for a greater value of the infec- 
tion probability than the corresponding value for the SIS 
model. This is expected because, in the SIS model, an S 
individual can be infected multiple times. 
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the square lattice occurs at cq = 0.1765005(10), consis- 
tent with, but much more precise than, previous work 
|19) . This critical value of c compares with c = 1/3 
for the SIR model on the 4-coordinated Bethe lattice, 
0.37753 for the contact process (or SIS model) on the 
square lattice, and 0.407254 for simultaneous infection of 
all neighboring S sites on the square lattice (site percola- 
tion). Having an accurate value of cq is useful for other 
studies of this critical state, such as studying the scaling 
of the average cluster size with L for finite systems [50] . 

Note added in revision. While this paper was under 
review, a paper appeared online j51j in which heterogene- 
ity of the transmissibility in a square-lattice SIR model 
was also considered. In that paper, the heterogeneity 
is formed by having classes of individuals each with dif- 
ferent fixed transmissibilities "01 corresponding to having 
different fixed infections times. In terms of the corre- 
lations we studied, each of these classes of individuals 
produce a binomial distribution for P^™^ given by Eq. 
( 10 1 with p replaced by V'i, and by taking a set of classes 



it is possible for example to reproduce the exponential in- 
fectivity distribution exactly. Our critical transmissibil- 
ity 0.53841 and its mean square deviation = 0.11575 
are consistent with the approximate criticality correla- 
tion these authors develop. More details will be given in 
a future publication. It should also be mentioned that 
that paper contains many references to epidemiological 
systems where a lattice-based SIR model is relevant. 
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IX. APPENDIX 

Consider that there are m nearest-neighbor S sites to 
a given I site, and four nearest neighbors total. Then we 
define: 

• ^nk ~ probability that exactly k distinct S 
sites are visited after n trials on the I site, k = 
0, . . . TO, for n = 0, 1, 2, . . . 00. 



We have provided further evidence that the asyn- 
chronous SIR model is in the universality class of stan- 
dard percolation through the behavior of the cluster size 
distribution. We showed that local correlations differ 
from those of standard percolation. By extensive numer- 
ical simulation of the cluster size distribution, we have 
shown that the transition in the SIR model defined on 



with Pg^"^ = 1; Pq^™^ = 0, = l,...,m. Given p{n) 



defined in ( 11 1, the net probability that k sites of type S 



are visited is given by 



(14) 



n=k 
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We can write recursion relations to find the ■ For 
example, for m = 1 we have 



p(i) _ 3 (1) 



p(l) _ 1 plJ-J I pV 



1 



,(1) 



(15) 
(16) 



The first expression states that the probability that the S 
site was not visited at the n-th trial is 3/4 the probability 
it was not visited in the previous trial, as three out of 
four of the neighbors are not an S site. Likewise, for the 
second expression, if the S was not visited in the previous 
trials, then in one out of four times the S will be chosen 
for the first time in the n-th trial, while if the S site was 
visited in some previous trial, then it will surely still have 
been visited in this trial. 

To solve the recursion relations, we use the generating 
functions, defined for general m by 



X 



(17) 



A straightforward application to ( TSpG ) yields 



1 - (3/4)x 



4(1 - a;)(l - (3/4)a;) 



(18) 
(19) 



According to (14), we have simply 



P^(™)=cGi"^(l-c) , 



(20) 



and we thus find 



p(i; _ 

J n — 



4c 



(1) 



1 + 3c 
1-c 
1 + 3c ' 



(21) 
(22) 



Likewise, for m = 2, we have 



(2) 



n,0 



(2) 



(2) 



n,l 



-P 



ri-1,0 
(2) 



3p(2) 

4 " 



p 



(2) 



n,2 



-P 



(2) 



ri-1.1 



P 



(2) 



n-1,2 



(23) 
(24) 
(25) 



and similarly one can write the recursions for m = 3 
and 4. In fact one can summarize them for all m by the 
general formulas 



P 



(m) 



4 — m 



P 



n,0 

(m) 



P 



(m) 



(26) 



TO + 1 



-p: 



(m) 



n,k ~ ^ ^ n-l,fc-l ^ ^ ^n-l,k i'^'^) 

for /c = 1, . . . , TO, and derive a general recursion relation 
for the averaged quantities for all to = 1, . . . , 4: 



P 



(m) 



4c 



TO + (4 — m)c 



(28) 



Pr ^ Plrl fc = l,...X29) 

TO - fc + (4 + A; - to)c ''"^ ' ' / 

Using these formulas, we obtain the numerical results 
given in Table [l] For to = 4 (the last row in that table), 
the explicit formulas are: 



P 



(4) 



(4) _ 4(1 - c)c 



p(^) = 



p 



(4) 



3 + c 

12(1 -c)2c 
(3 + c)(2 + 2c) 

24(1 - cfc 
(3 + c)(2 + 2c)(l + 3c) 

6(1 -c)^ 
(3 + c)(2 + 2c)(l + 3c) 



(30) 
(31) 

(32) 
(33) 
(34) 



Note that it is possible to write explicit formulas for 

all the P^™^ , such as 

(3) _ 4" ~ 3(3")+ 3(2") -1 
where the coefficients in some cases are related to Stirling 



numbers of the second kind. Also, we can generalize ( 28 



29 1 to all coordination numbers z simply by replacing all 



4's in those equations by z's. Then one can show that 



P 



(m) 



P 



(„) ^ n{l - c) („) 



which satisfies P^^^ + P^"' = 1. 



for n = 1,...,TO, and taking n — m, one can verify 
directly that Er=o ^i"^ = 1- 
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